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1.  INTRODUCTION 


The  one-dimensional  Navy  Standard  Surf  Model  (NSSM,  or  SURF  3.2)  has  been 
shown  to  be  very  robust  (Hsu,  et  ah,  2002),  but  it  can  produce  inaccurate  wave  and 
longshore  current  estimations  for  areas  with  complicated  bathymetry.  Since  SURF 
assumes  parallel  bottom  contours  in  the  surf  zone,  it  cannot  account  for  longshore 
variations  of  bathymetry  or  forcing.  Only  2D  or  quasi-3D  nearshore  models  should  be 
used  for  such  cases.  The  Delft3D  modeling  system,  developed  by  Delft  Hydraulics 
(Roevink  and  Banning,  1994),  is  a  complete  coastal  hydrodynamic  modeling  system, 
capable  of  simulating  hydrodynamic  processes  due  to  waves,  tides,  rivers,  winds  and 
coastal  currents.  It  uses  finite  difference  numerical  techniques  and  can  be  run  in 
Cartesian  or  spherical  coordinates  and  for  regular  or  curvilinear  grids.  Delft3D  has  been 
evaluated  and  validated  by  Elias  et  al.  (2000)  and  Morris  (2001).  Recently,  Hsu  et  al. 
(2006)  validated  the  model  using  Duck94  and  Santa  Barbara  data.  Additional  validation 
using  SandyDuck  (Duck,NC)  and  SURF04  (Sardinia/Italy)  data  sets  are  being  conducted. 

Delft3D  produces  two-dimensional  time  dependent  forecasting  output  for  many 
nearshore  wave  and  flow  parameters.  But  it  does  not  produce  the  operational  surf 
forecasting  parameters  as  specified  in  the  Joint  Surf  Manual.  Using  Delft3D  output,  this 
report  describes  the  software  which  generates  the  standard  surf  parameters  including 
maximum  and  significant  breaker  height,  breaker  type  statistics,  percent  of  breaking,  surf 
zone  width,  number  of  surf  lines  and  modified  surf  index  (MSI).  Subroutines  from 
SURF  3.2  are  adapted  to  compute  these  parameters. 

2.  DeIft3D  OUTPUT  FILES 

As  described  in  the  Delft3D-FLOW  manual  (Delft  Hydraulics,  2005),  the  results  of  a 
Delft3D-FLOW  computation  are  stored  in  four  types  of  files: 

•  communication  file:  <com-mnid.def>  and  <com-mnid.dat>. 

•  history  file:  <trih-runid.def>  and  <trih-mnid.dat>. 

•  map  file:  trim-mnid.def>  and  <trim-mnid.dat>. 

•  drogue  file:  <trid-runid.def>  and  <trid-mnid.dat>. 

These  files  use  the  binary  NEFIS  format,  therefore  is  not  easily  readable.  MATLAB  routines 
provided  by  Delft3D  are  used  in  this  software  development  to  read  the  data.  All  necessary  input 
(cross-shore  data  of  location,  depth,  significant  wave  height,  wave  angle  and  longshore  current)  to 
derive  the  standard  surf  parameters  can  be  obtained  from  the  “com”  (communication)  file.  If  the 
roller  option  is  tuned  on,  then  significant  wave  height  needs  to  be  derived  from  the  wave  energy 
term  in  the  “trim”  (map)  file.  When  roller  is  on,  radiation  stresses  and  gradients  of  these  stresses 
are  computed  in  Delft3D  based  on  the  wave  energy  and  roller  energy  replacing  the  conventional 
wave  forces  as  derived  from  the  WAVE  (SWAN)  model.  Only  the  wave  direction  from  an  initial 
SWAN  mn  is  still  used.  It  should  be  noted  that  the  currents  derived  from  trim  file  represent  the 
Eulerian-frame  (fixed  location)  velocity  whereas  those  from  com  file  represent  Lagrangian-frame 
(following  tagged  water  particles)  velocity.  But  as  far  as  longshore  current  is  concerned,  currents 
derived  from  both  files  are  the  same. 
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3.  SOFTWARE  DESCRIPTION 


To  help  describe  the  objective  of  the  software,  a  sample  surf  summary  for  a  test 
case  from  SURF  3.2  is  presented  in  Appendix  A.  The  first  part  of  the  summary  is  the 
input  information,  and  the  part  starting  with  “Coded  Surf  Forecast  Follows”  is  what  this 
software  will  produce  from  DelftSD  results.  Fig.  1  shows  the  flow  chart  of  the  MATLAB 
code  to  convert  DelftSD  binary  output  to  ASCII  format. 


Fig.  1  -  Flow  chart  for  Matlab  delft3d2ascii.m 


Two  types  of  output  files  are  generated.  The  first  one  is  called  dSd.list  which 
shows  the  list  of  data  filenames  at  different  alongshore  locations.  A  partial  sample  is 
listed  here: 

d3d_0.dat 

d3d_30.dat 

d3d_60.dat 

d3d_90.dat 


d3d_1710.dat 

d3d_1740.dat 

The  number  in  the  filenames  represents  the  longshore  positions  of  the  cross-shore  data, 
e.g.,  d3d_90.dat  is  for  data  at  90-m  location  relative  to  the  model  grid  origin.  The  second 
file  consists  of  cross-shore  data  at  a  particular  alongshore  location  and  is  partially  listed 
below: 
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X 

900.00 

890.00 


depth 

8.26 

8.19 


wave  longshore  wave 
height  current  angle 


2.34  -0.07  -7.8 
2.34  -0.07  -7.8 


310.00 

3.83 

1.92 

-0.50 

-3.1 

300.00 

3.60 

1.88 

-0.51 

-2.7 

290.00 

3.32 

1.83 

-0.51 

-2.3 

280.00 

2.95 

1.76 

-0.49 

-1.7 

270.00 

2.50 

1.64 

-0.46 

-1.0 

Fig.  2  shows  the  flow  chart  of  the  main  program  where  the  subroutine  names  are 
listed  in  parenthesis.  All  subroutines  are  adopted  from  SURF  3.2.  In  additional  to  these 
two  types  of  files,  another  input  file  called  input. par  consists  of  beach  orientation,  peak 
wave  period  offshore,  wind  direction,  wind  speed  and  output  interval.  Similar  to  SURF 
3.2,  beach  orientation  angle  is  defined  as  the  compass  heading  towards  beach.  The 
output  interval  represents  the  alongshore  distance  between  adjacent  output.  It  is  noted 
that  all  input  files  are  in  MKS  unit,  i.e.  in  meter  etc. 


Fig.  2  -  Flow  chart  for  main  Fortran  program 
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Two  types  of  output  files  are  generated  by  this  software.  The  first  type  of  files  is 
named  with  an  “out”  extension: 

d3d_0.out 

d3d_30.out 

d3d_60.out 

d3d_90.out 


They  correspond  to  the  surf  summary  for  alongshore  locations  from  those  input  files  as 
listed  in  d3d.list,  e.g.  d3d_540.dat  is  for  summary  for  input  d3d_540.dat.  Identical  in 
content  to  the  SURF  3.2  summary  as  listed  in  Appendix  A,  a  sample  “out”  file  from 
Delft3D  results  is  listed  here: 


*****  *****  Coded  Surf  Forecast  Follows  ***** 


Significant  Breaker  Height 
Maximum  Breaker  Height 
Dominant  Breaker  Period 
Dominant  Breaker  Type 
(  98%  Spilling,  2%  Plunging,  0% 
Breaker  Angle  (toward  left  flank) 
Littoral  Current  (toward  left  flank) 
Number  of  Surf  Lines 
Surf  Zone  Width 
Average  Wave  Length 
Wind  Speed 
Wind  Direction 


alfa  = 
bravo  = 
Charlie  = 
delta  = 

Surging) 

echo  = 
foxtrot  = 
golfl  = 
golf2  = 

hotell  = 
hotel2  = 


Modified  Surf  Index  =  7.2 


★  ★  ★  ★  ★ 


4.9  ft 
7.5  ft 
8 . 0  sec 
Spilling  Surf 


12.9  deg 
1.7  kts 
4.2 

538.0  ft 
126.7  ft 
10.0  kts 
46.0  deg 


Units  are  in  feet  as  in  the  SURF  3.2.  The  second  type  of  output  file  is  called  msi.dat 
consisting  of  MSI,  surf  zone  width  and  it’s  alongshore  distance.  Many  other  surf 
parameters  can  easily  be  added  on  the  output  list. 
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4.  SAMPLE  RESULTS 


A  simulation  of  DelftSD  was  conducted  for  a  Duck94  case  at  1600  EST,  10 
October  1994.  The  DelftSD  output  was  processed  by  the  new  software.  Fig.  3  shows  the 
computed  MSI  (from  the  output  file:  msi.dat)  as  a  function  of  longshore  locations.  The 
straight  line  at  MSI  =  8  represents  operational  limit  of  landing  craft  mechanized  (LCM-6, 
LCM-8)  and  landing  craft,  air  cushion  (LCAC). 


Fig.  3  -  Sample  plot  of  MSI  for  a  Duck94  test  case. 

In  Fig.  4,  the  edge  of  surf  zone,  derived  from  the  computed  surf  zone  width,  is  plotted 
with  longshore  current  over  depth.  The  white  line  corresponds  to  locations  where  10 
percent  of  waves  are  breaking.  The  yellow  arrow-head  represents  1  knot  (0.51  m/s)  of 
current. 
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Surf  Zone  and  Current  over  Depth,  Duck/NC,  Input  Wave  2  m 


Fig.  4  -  Sample  plot  showing  the  edge  of  surf  zone  with  longshore  eurrents  over  depth. 

5.  FUTURE  PLANS 

The  present  software  applies  to  Cartesian  grids  where  the  coastline  is  almost 
straight,  so  that  beach  angle  remains  about  the  same  for  the  entire  modeling  region. 
Future  refinement  of  the  software  should  extend  its  application  to  curvilinear  grids.  The 
present  operational  Perl  script,  for  mnning  DelftSD  using  DIOPS  (Allard,  et  ah,  2005) 
output,  should  also  incorporate  the  present  software,  so  that  standard  Navy  surf  output 
can  be  generated  automatically. 
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APPENDICES 


Appendix  A  -  Sample  SURF  3.2  output  showing  the  standard  Navy  surf  output 
Test  Case  1  Output  File  -  case  Lout 

*****  *****  Surf  Forecast  *****  ***** 


Navy  Standard  Surf  Model  Version  3.2 
Date  and  Time  of  Forecast:  01/01/2001  0100 

Session  Logged  to  file  easel. out 
Landing  Zone  Name  =  case  7  beach 

Sight  Line  =  0.0  deg 

Equilibrium  Beach  Sediment  =  medium  sand 
Wave  Input  Depth  =  25.0  ft 

Sea  Height,  Period,  Direction  =  0.0  ft,  0.0  sec. 

Swell  Height,  Period,  Direction  =  3.0  ft,  10.0  sec. 

Wind  Speed  =  10.0  kts 

Wind  Direction  =  240.0  deg 

Tide  Level  =  -1.0  ft 


0 . 0 
160 . 0 


deg 

deg 


Internally  Generated  Spectrum  Used 


starting  Depth  = 

24  . 

0 

ft 

Output  Interval  = 

5. 

0 

ft 

Computational  grid  spacing  = 

2  . 

0 

ft 

Significant  Wave  Height  Offs 

hore 

= 

3. 

,  0 

ft 

Peak  Period 

= 

10  . 

,  0 

sec 

Average  wave  direction 

= 

-20  . 

,  0 

deg 

Percent  of  Breaking  Waves  is 

less 

than 

5. 

,  0 

"0 

starting  depth. 


■k  -k  -k  -k  -k  -k  -k  -k  -k  -k 


alf  a 
bravo 
Charlie 
delta 

Surging) 


Coded  Surf  Forecast  Follows 
Significant  Breaker  Height 
Maximum  Breaker  Height 
Dominant  Breaker  Period 
Dominant  Breaker  Type 
(  84%  Spilling,  16%  Plunging,  0% 

Breaker  Angle  (toward  left  flank) 

Littoral  Current  (toward  left  flank) 

Number  of  Surf  Lines 
Surf  Zone  Width 
Average  Wave  Length 
Wind  Speed 
Wind  Direction 


■k  -k  -k  -k  -k  -k  -k  -k  -k  -k 


3.6  ft 
5.5  ft 

=  10.0  sec 

=  Spilling  Surf 


echo  = 

11 . 

9 

deg 

foxtrot  = 

1 . 

5 

kts 

golfl  = 

2  . 

0 

golf2  = 

246. 

0 

ft 

= 

120  . 

9 

ft 

hotell  = 

10  . 

0 

kts 

hotel2  = 

240  . 

0 

deg 

Modified  Surf  Index  = 


7.2 


*****  *****  Detailed  Surf  Output  Follows  *****  ***** 


Indx 

Dist 

Water 

Sig  Brkr 

Max  Brkr 

Prent 

Brkr 

Littoral 

Offshore 

Depth 

Height 

Height 

Brkng 

Angle 

Current 

(ft) 

(ft) 

(ft) 

(ft) 

waves 

(deg) 

(kts ) 

1 

302 . 4 

9.3 

3.57 

5.46 

5.2 

-12 . 7 

0 . 12 

2 

297.4 

9.1 

3.57 

5.47 

5.5 

-12 . 6 

0 . 07 

3 

292 . 4 

9.0 

3.57 

5.47 

5.9 

-12 . 5 

0 .01 
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Appendix  B  -  MATLAB  delft3d2ascii.m  to  read  DelftSD  NEFIS  format  output  files 


%function  []  =  delft3d2ascii  (runid, datetime, ctitle, roller ) 

%  Read  the  DelftSD  com  and  trim  flies  to  produce  plots  and  ascii  output 
flies 

%  Developed  by  E.  Rogers,  J.  Dykes  and  L.  Hsu,  Code  7322,  NRL 
%  The  following  parameters  are  placed  here  to  be  talloured  to  your 
%  particular  problem  and  preferences. 

path (path, ' /home/utah/hsu/delft3d/2 . 11 . 00/d3dmatlab ' ) 
runld= ' du8 ' ; 
datetlme=' 1994101016 ' 
ctltle='Duck94  '  ; 

roller  =1;  %when  roller  Is  used  In  *.mdf  file,  otherwise  roller=0 
%  When  roller  Is  used,  wave  height  should  be  derived  from  trim  Instead 
com  file 

thln_s  =4;  %  plot  thinning  of  arrows  In  current  plot 
thln_w  =4;  %  plot  thinning  In  arrows  In  wave  plot 

quality  =  '75';  %  jpeg  quality  level,  affects  figure  size  and  quality 

hscale  =  75;  %wave  direction  arrow  scale 

vscale  =  125;  %veloclty  arrow  scale 

vlegend  =  0.5;  %arrowhead  legend  velocity 

range  =  [-18:3:0];  %depth  contour  Interval  and  range 

sc  =  0.25;  %  contour  smoothing,  fine  grid  uses  a  small  number,  coarse 
grid  uses  a  large  number 
fz  =  12;  %font  size 

%  For  cropping 

%xstart  =  50;  %  staring  Index 

%nxend  =  171;  %  ending  Index 

%  Com  file  Is  opened  and  relevant  Info  Is  extracted, 
filename  =  sprlntf  ('%s%s',  'com-',  runid,  ' . dat ' ) 

Nfs  =  vs_use  (f llename, ' quiet ') ; 

fllename2  =  sprlntf  ('%s%s',  'trim-',  runid,  '.dat') 

Nfs2  =  vs_use  (fllename2, ' quiet ') ; 

%vs_dlsp  (Nfs);  %  show  all  stored  parameters 
%ntcurs  =  vs_get  (Nf s ,' KENMNT ',' NTCUR ',' quiet ') ; 
ref_date  =  vs_get  (Nfs,  ' PARAMS ' ,  'ITOl',  'quiet'); 
ref_tlme  =  vs_get  (Nfs,  'PARAMS',  'IT02',  'quiet'); 
tscale  =  vs_get  (Nfs,  'PARAMS',  'TSCALE',  'quiet'); 
tlmcur  =  vs_get  (Nf s ,' CURTIM ',' TIMCUR ',' quiet ') ; 

If  Iscell  (tlmcur)  ==  0 
ntlme  =  tlmcur; 
else 

ntlme  =  cell2mat  (tlmcur) ; 
end 

tau  =  tscale  *  ntlme ( 1 ) /3600 ; 
cref_date  =  num2str  (ref_date); 
cref_tlme  =  num2str  (ref_tlme); 
yyyy  =  str2num  (cref_date  (1:4)); 
mm  =  str2num  (cref_date  (5:6)); 
dd  =  str2num  (cref_date  (7:8)); 
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hh  =  str2nuin  (cref_tiine)  ; 

ntime  =  datenum  (yyyy,  mm,  dd,  hh,  0,  0)  +  tau/24; 

[vyyyy,  vmm,  vdd,  vhh]  =  datevec  (ntime); 

valid_date  =  sprintf  ( ' %d%0 . 2d%0 . 2d '  ,  vyyyy,  vmm,  vdd); 

%CLH  valid_hour  =  sprintf  ('%0.2d',  vhh); 

%somehow  getting  date  and  time  from  above  is  not  always  correct,  so 

datetime  (an  input  from  script)  instead 

valid_date=datetime (1:8)  ; 

vhh  =str2num (datetime ( 9 : 10 ))  ; 

yymmddhh=dat  et ime (3:10) ; 

valid_hour  =  sprintf  ('%0.2d',  vhh); 

%  Get  currents 

%data  =  vs_get  (Nf s ,' CURTIM U1 quiet ') ; 
data  =  vs_get  (Nf s2 ,' map-series U1 quiet ') ; 
nt  =  size  (data,  1);  %  nt  is  for  getting  the  last  data  set 
%  if  only  one  data  set  (one  time  step)  used,  the  data  array  is  no 
longer  a  cell  array, 
if  iscell  (data)  ==  0 
USPD  =  data ' ; 
else 

USPD_CELL  =  data(nt); 

USPD  =  USPD_CELL{ 1, 1 } ' ; 
end 

data  =  vs_get  (Nfs2, 'map-series ', 'VI quiet ') ; 

%data  =  vs_get  (Nf s ,' CURTIM ',' VI ',' quiet ') ; 
if  iscell (data)  ==  0 
VSPD  =  data ' ; 
else 

VSPD_CELL  =  data(nt); 

VSPD  =  VSPD_CELL{ 1, 1 } ' ; 
end 

%  Get  wave  parameters 

%  when  using  roller,  get  wave  energy,  then  convert  to  wave  height 
if  roller  ==  1 

data  =  vs_get  (Nfs2,  ' map-rol-series ' ,  'EWAVEl',  'quiet'); 
if  iscell  (data)  ==  0 
weng  =  data ' ; 
else 

weng_cell  =  data(nt); 
weng  =  weng_cell { 1 ,  1 }  '  ; 
end 

%  Change  engery  to  wave  height 
whrms=sqrt (8*weng/ 9.81/1025)  ; 

else 

whrms  =  vs_get  (Nfs,  'WAVTIM',  ' HRMS ' ,  'quiet') '; 
end  ;  %  end  of  roller  logic 

wdir  =  vs_get  (Nfs,  'WAVTIM',  'DIR',  'quiet') '; 

DPO  =  vs_get  (Nfs,  'INITBOT',  'DPO',  'quiet')'; 

X2D  =  vs_get  (Nfs,  'GRID',  'XCOR',  'quiet')'; 


use 
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Y2D  =  vs_get  (Nfs,  'GRID',  'YCOR',  'quiet')'; 

%  Get  rid  of  extra  row  and  column  put  in  by  DelftSD 
nx  =  size  (X2D,  1)  -  1; 
ny  =  size  (X2D,  2)  -  1; 

DPO  =  -DP0(l:nx,  l:ny);  %  set  depth  as  negative 
X2D  =  X2D (1 :nx,  1 :ny) ; 

Y2D  =  Y2D (1 :nx,  1 :ny) ; 
x=X2D ( : , 1) ; 
y=Y2D  (1,  :  )  ; 

USPD=USPD (1 :nx, 1  :ny)  ; 

VSPD=VSPD ( 1 : nx, 1 : ny ) ; 
whsig=sqrt (2 ) *whrms ( 1 : nx, 1 : ny) ; 
wdir=wdir ( 1 : nx, 1 : ny ) ; 

%  Crop  boundaries  from  m  =  nxstart  to  nxend 
%CLH  no  cropping  nx  =  nxend  -  nxstart  +  1; 

%Y2D  =  Y2D (nxstart : nxend,  :); 

%X2D  =  X2D (nxstart : nxend,  :); 

%USPD  =  USPD (nxstart : nxend,  :); 

%VSPD  =  VSPD (nxstart : nxend,  :); 

%DP0  =  DPO (nxstart : nxend,  :); 

%whrms  =  whrms (nxstart : nxend,  :); 

%wdir  =  wdir (nxstart : nxend,  :); 

%  thinning 
it  =  0 ; 

%skip  near  the  boundary 
jskip  =  0; 

for  i  =  l:thin_s:nx 
it  =  it  +  1; 
jt  =  0; 

for  j  =  (jskip+1) :thin_s: (ny  -  jskip) 
jt  =  jt  +  1; 

Y2D_s (it,  jt)  =  Y2D(i,  j); 

X2D_s (it,  jt)  =  X2D(i,  j); 

USPD_s(it,  jt)  =  USPD(i,  j); 

VSPD_s(it,  jt)  =VSPD(i,  j); 
end 
end 

it  =  0 ; 

for  i  =  l:thin_w:nx 
it  =  it  +  1; 
jt  =  0; 

%  for  j  =  (jskip+1) :2*thin_w: (ny  -  jskip) 
for  j  =  (jskip+1) :thin_w: (ny  -  jskip) 
jt  =  jt  +  1; 

Y2D_w(it,  jt)  =  Y2D(i,  j); 

X2D_w(it,  jt)  =  X2D(i,  j); 

whx(it,  jt)  =  whrms (i,  j)  *  cos  (wdir(i,  j)  *  pi/180); 

why(it,  jt)  =  whrms (i,  j)  *  sin  (wdir(i,  j)  *  pi/180); 

end 

end 

figure  (1),  clf,  hold  off 
axis  equal 
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h  =  pcolor  (x,  y,  DPO  '  ) ; 
axis  equal 

set  (h,  'linestyle',  'none') 
shading  interp;  %  smooth  looking 
%axis  ([min(x)  max(x)  min (y)  max(y)]) 
newmap= jet ; 

% light  jet =newmap (12:58,  :) ; 

colormap ( jet ) ; 

caxis  ( [min (min (DPO ) )  0]) 

colorbar 

xlabel  ('x  (m) ', ' f ontsize ' , f z ) 

ylabel  ('y  (m) ',' f ontsize ', fz ) 

titlel  =  sprintf  (  ' %s%s%s%s%s%s%s  '  ,  'Current  over  Depth,  ',  ctitle, 

',  valid_date,  '  ',  valid_hour,  'EST'); 

title  ([titlel]  ,' f ontsize ', 10 ) 

set  (gca, ' f ontsize ', fz ) 

hold  on 

quiver  (X2D_s,  Y2D_s,  USPD_s*vscale,  VSPD_s*vscale,  0,  'k-'); 

set  (findobj  ('type',  'line'),  'linewidth',  0.03) 
hold  on 

xd  =  (max (x) -min (x) ) /10+min (x) ; 
yd  =  (max  (y ) -min  (y )  )  / 1 6+min  (y )  ; 
xd2  =  (max (x) -min (x) ) /11+min (x) ; 
yd2  =  (max  (y ) -min  (y )  )  / 1 1+min  (y )  ; 

%  current  legend 

quiver (xd, yd, vscale*vlegend, 0 . , 0 , ' m- ' ) 

text (xd2, yd2, '0.5  m/ s ' , ' font size ',11, 'color', 'm') 

figure  (2) 
axis  equal 

h  =  pcolor  (x,  y,  whsig'); 
axis  equal 

set  (h, 'linestyle', 'none') 

shading  interp;  %  for  smooth  looking 

colormap ( jet ) 

colorbar 

caxis  (  [0  max (max (whsig) ) ] ) 
xlabel  ('x  (m) ',  'fontsize',  fz) 
ylabel  ('y  (m) ',  'fontsize',  fz) 

title2  =  sprintf  ( ' %s%s%s%s%s%s%s ' ,  'Sig  Wav  Hgt  and  Depth  Cont,  ', 

ctitle,  '  ',  valid_date,  '  ',  valid_hour,  'Z'); 

title  ([title2]  ,' fontsize ', 10 ) 

set  (gca,  'fontsize',  fz) 

hold  on 

quiver  (X2D_w, Y2D_w, whx*hscale, why*hscale, 0, ' k- ' ) ;  %  wave  vector 
set  (findobj ( ' type ' , ' line ' ) , ' linewidth ',0.03) 
hold  on; 

xf  =  linspace  (min(x),  max(x),  nx*sc)  ; 
yf  =  linspace  (min (y) ,  max (y) ,  ny*sc) ; 

[X,Y]  =  meshgrid  (x,y); 

[XI, YI]  =  meshgrid  (xf,yf); 

ZI  =  interp2  (X,  Y,  DPO',  XI,  YI,  'cubic'); 

[C,h]  =  contour  (xf,  yf,  ZI,  range,  'm — '); 
set  (h,  'linewidth',  0.03) 
hh=c label (C, h, ' Label Spacing ',216) ; 
set  (hh,  'color',  'm' ,' fontsize ' ,  11) 
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vduc  '  yyminddhh]  ; 


str  =  ['print  -fl  -dpng  ' 
disp  (str) 
eval  (str) 

%str  =  ['print  -f2  -dpng  '  'w'  datetime  ]; 

str  =  ['print  -f2  -dpng  '  'wduc'  yymmddhh  ]; 
disp  (str) ; 

%  If  the  coordinates  were  in  UTM,  this  tabular  output  would  be  suitable 
for  ingest  into  ArcView. 

%fileout  =  sprintf  ('%s%s%s%s',  ctitle,  'table',  ' . dat ' ) ; 

%fout  =  fopen  (fileout,  'w'); 

%fprintf  (fout,  '%9s  %9s  %9s  %9s  %9s  %9s  %9s\n',  'X',  'Y',  'Depth(in)  ', 

'Hsig(in)',  'Dir(Cart)',  'U',  'V'); 

%in  =  0; 

%f or  iy  =  1 : ny 
%  for  ix  =  1 : nx 
%  if  (DP0(ix,iy)  <  0) 

%  fprintf  (fout,  '%9.1f  %9.1f  %9.1f  %9.1f  %9.1f  %9.4f  %9.4f\n', 

X2D(ix,iy),  Y2D(ix,iy),  DP0(ix,iy),  whs ig ( ix, iy ) ,  wdir(ix,iy), 
USPD(ix,iy),  VSPD ( ix, iy )  )  ; 

%  else 

%  fprintf  (fout,  '%9.1f  %9.1f  %9.1f  %9.1f  %9.1f  %9.4f  %9.4f\n', 

X2D(ix,iy),  Y2D(ix,iy),  DP0(ix,iy),  NaN,  NaN,  NaN,  NaN) ; 

%  end 

%  end 
%end 

%f close  (fout ) ; 

f idl=f open ('dSd.list',  'w')  ; 

f or  j  j  =  1  :  2 :  ny 
ydist  =  15* ( j  j-1 )  ; 

%  change  cartesian  to  surf  wave  angle  convention 
tmpwdir ( : , j  j ) =- (wdir (:,  jj)-180)  ; 
for  ii=l : nx  ; 

%  fix  for  wave  angle  beyond  computation,  i.e.  on  land 
if  (tmpwdir (ii, j j )  ==  180); 
tmpwdir (ii, jj)=0.0; 
end 
end 

%only  save  the  data  when  depth  is  larger  than  0,  i.e.  wet 
ztmp=-DP0 ( : , j  j ) ; 

L=find(ztmp  >  0); 

output=  [ ( f lipud (X2D (L, j j ) ) ) ' ;  (flipud(- 
DPO  (L, j  j ) ) )  ' ;  (f lipud (whsig (L,  j  j ) ) )  '  ; 

( f lipud (VSPD(L,jj)))  ';  ( f lipud (tmpwdir (L, j  j ) ) )  ' ] ; 
prof ileout=sprintf  ('%s',  ' d3d_ ' , num2str (ydist ) ,  '.dat'); 

f id=f open (prof ileout ,  ' w ' )  ; 

fprintf (fid, ' %8 . 2f  %8.2f  %8.2f  %8.2f  %8.1f\n  ', output); 
fclose (fid) ; 

fprintf(fidl,  '%s\n  ',  prof ileout ) ; 
end 

fclose (fidl) ; 
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Appendix  C  -  Main  Fortran  Program  surfdelftSd.f 
(All  subroutines  are  adapted  from  SURF  3.2  and  are  not  listed  here) 


C  This  program  is  developed  to  produce  standard  Navy  surf  parameters 
C  from  DelftSD 

C  All  subroutines  are  adapted  from  SURF  3.2. 

C  Programmed  by  Larry  Hsu  and  validated  by  Ted  Mettlach 
C  For  questions,  contact  Larry  Hsu,  code  7322,  NRL, 
hsu@nrlssc. navy . mil 

C  Input  files  are  produced  from  delf t2ascii . m  (MATLAB) :  d3d.list  and 
C  d3d_ydist.dat  files.  Additional  input  file  is  input. par  which 
provides 

C  beach  orientation,  peak  wave  period  offshore,  wind  direction,  wind 
speed 

C  and  output  interval. 

C  Output  files  are  msi.dat  and  d3d_ydist . out  files  which  give  surf 
summary . 

c  local  variables  from  SURF  3.2: 
c  alfa  real  significant  breaker  height 

c  bravo  real  maximum  breaker  height 

c  chrlie  real  dominant  breaker  period 

c  dangle  real  angle  between  directional  bins 

c  depname  char*40  depth  profile  file  name 

c  dsea  real  input  direction  for  sea  contribution 

c  dstart  real  input  starting  depth 

c  dswell  real  input  swell  direction  for  internally 

c  generated  spectrum 

c  dxyl (points)  real  corresponding  depths  with  no  tide 

c  echo  real  breaker  angle 

c  ehsig  real  significant  wave  height  from  directional 
c  spectrum 

c  esowm  (dirnum, freqnum)  real  directional  wave  spectrum 

c  file_dat  char*40  output  file  name  * . dat 

c  file_in  char*40  input  filename 

c  file_out  char*40  output  file  name  * .  out 

c  file_tmp  char*40  temporary  file 

c  file_spc  char*40  output  file  name  for  refracted  wave  spectrum 

c  foxtrt  real  longshore  current  speed  and  direction 

c  fracname  char*40  wave  refraction  file  name 

c  freq  (freqnum)  real  input  wave  spectrum  center  frequencies 

c  freq_rs ( freqnum)  real  frequencies  assoc'd  w/  input  refracttion  & 
c  shoal  matrices 

c  freql  (freqnum)  real  beginning  frequency  bin  values 

c  freq2  (freqnum)  real  ending  frequency  bin  values 

c  gamma2  real  beach  orientation,  compass  leading  directly  toward  beach 
c  golfl  real  number  of  surf  lines 

c  golf2  real  surf  zone  width 

c  hsea  real  input  significant  wave  height  for  sea 

c  contribution  to  pierson  moskowi 

c  hswell  real  input  significant  wave  height  for 

c  internally  generated  spectrum 

c  iday  integer  input  day 
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c  idirec  integer  number  of  direction  bins  in  the  input  spectrum 
c  idirec_rs  integer  number  of  directions  assoc 'd  w/  input 
c  refraction  &  shoal  matrices 

c  ifreq  integer  number  of  frequency  bands  in  the  input  spectrum 
c  ifreq_rs  integer  no.  of  freqs  assoc 'd  w/  input  refraction 

c  &  shoal  matrices 

c  igamma  integer  beach  orientation  rotated  90°  from 
c  original  heading  toward  beach 

c  ihour  integer  input  hour 

c  ihtll  real  wind  speed  coded  surf  forecast  value 
c  ihtl2  real  wind  direction 
c  imin  integer  input  minute 
c  imonth  integer  input  month 
c  iyear  integer  input  year 

c  jgamma  integer  temporary  value  set  to  beach  orientation 
c  line  char*80  temporary  variable  used  to  read  lines 
c  lin_stress  logical  longshore  current  solution  (true  or  false) 
c  Indname  char*40  input  landing  zone  name 
c  nnn  integer  number  of  points  in  the  input  depth  array 
c  pet  (4)  real  percent  of  different  breaker  types 


c 

pet 

(1) 

=  spilling 

c 

pet 

(2) 

=  plunging 

c 

pet 

(3) 

=  surging 

c 

pet 

(4) 

=  total 

c  period ( freqnum)  real  period  array  (1/frequency) 

c  psea  real  input  wave  period  for  sea  contribution  to 
c  pierson  spectrum 

c  pswell  real  input  swell  period  for  internally  generated 
c  spectrum 

c  roller  logical  roller  usage  (true  or  false) 
c  self_st  char*l  self  start  flag  (yes  or  no) 
c  slope  real  bottom  slope 

c  spectra  logical  does  input  spectra  exist?  (true  or  false) 
c  spefile  char*40  selected  wave  spectrum  file  name 
c  spedepth  real  depth  at  input  wave  spectrum 

c  surfy  logical  significant  wave  heights  greater  than 
c  0.5  ft?  (true  or  false) 

c  tide  real  input  tide  level 

c  wdir  real  input  wind  direction  compass  heading 

c  from  which  wind  comes 

c  wspd  real  input  wind  speed 

c  xcoeff (dirnum,  freqnum)  real  shoaling  coefficient  matrix  adjusted  to 
c  xfrom  and  freq 

c  xcoeff_rs (dirnum, freqnum)  real  shoaling  coefficient  matrix  input  from 
f racname 

c  xdelt  real  surf  zone  output  interval 
c  xdelt_gr  real  self-adjusting  cross-shore  grid  step 

c  xfrom (dirnum)  real  direction  array,  direction  from  which  wave  energy 
comes 

c  xfrom_rs  (dirnum)  real  direction  array  assoc 'd  w/  input 
c  refraction  &  shoal  martices 

c  xtheta (dirnum, freqnum)  real  refraction  angle  matrix  adjusted  to  xfrom 
and  freq 

c  xtheta_rs (dirnum, freqnum) real  refraction  angle  matrix  input  from 
f racname 

c  xxl (points)  real  adjusted  cross-shore  distances  from  depth  profile 
c  ydepth  char*l  input  depth  profile  used?  (yes  or  no) 
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c  ydetail  char*l  detailed  output?  (yes  or  no) 

c  yrefrac  char*l  is  refraction  considered  in  analysis?  (yes  or  no) 
c  ystr  char*l  is  strait  coast  refraction  used?  (yes  or  no) 
c 

program  surfdelftSd 
c 

implicit  none 
include  ' common. inc' 
c 

integer  ifreq,  idirec,  nnn,  igamma 
integer  ifreq_rs,  idirec_rs 

integer  iyear,  imonth,  iday,  ihour,  imin,  jgamma 
integer  spe_type 
c 

real  xxl (points),  dxyl (points ) ,  xf rom (dirnum) ,  pet (4) 
real  xfrom_rs (dirnum) 

real  f req ( f reqnum) ,  period ( freqnum) , 
freql (freqnum) , freq2 (freqnum) 
real  freq_rs (freqnum) 

real  esowm (dirnum,  freqnum),  xcoeff (dirnum,  freqnum) 
real  xcoeff_rs (dirnum,  freqnum) 
real  xtheta (dirnum,  freqnum) 
real  xtheta_rs (dirnum,  freqnum) 

real  alfa,  bravo,  chrlie,  echo,  foxtrt,  golfl,  golf2 
real  ehsig,  dangle,  slope,  gamma2,  dstart,  xdelt,  xdelt_gr 
real  hsea,  psea,  dsea,  hswell,  pswell,  dswell 
real  wspd,  wdir,  tide,  ihtll,  ihtl2,  spedepth 
c 

character*!  ydepth,  ydetail,  self_st,  yrefrac,  ystr 
character*40  file_in,  file_out,  fracname,  depname,  file_dat 
character*40  file_tmp,  spefile,  file_spc 
character*40  Indname 
c 

logical  roller,  lin_stress,  spectra,  surfy 

real  tmp (points),  duml,  dum2 

real  distmax 

real  fqd,  per,  dp 

real  dxy (points) 

real  xtemp (points ) ,  xktemp (points ) ,  htemp (points ) ,  ptemp (points ) 
real  thetatemp (points ) ,  thetamin,  thetamax 

real  ebtemp (points ) ,  along (points ) ,  blong (points ) ,  v (points) 
real  dong  (points )  ,  rk(points,  4),  bl  (points) 
real  suml,  width,  hlmax,  h2max,  vmin,  vmax 
real  df,  ftsq2msq,  temp 

integer  iimax,  wid_ii,  nnnl,  nnn2,  nnn3,  nnn4 
real  hrms,  xk 

integer  j_ii,  j,  k  ,  i  ,  j_ii2 
real  xshift 
character*40  file_in2 

real  tmph (points ) ,  tmptheta (points ) ,  tmpv (points) 
real  xxin (points ) ,  zzin  (points) 
real  xkd,  10,  xof f , p (4 ) , dum 
logical  surf,  brklO 

integer  ii,  ndepth,  fend,  irun,  iline,  iydist 
character*20  tmpfile(80),  cydist,  cdum 

real  tmpmsi(80),  tmpzone(80),  srfmod,  tmppet (80, 4) ,  tmpydist(80) 
real  xmin,  tmpxmin(80) 
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.false . / 


c  setting  model  options 

data  roller,  lin_stress  /.true., 
c 

ififi<-kififi<ifififi<-kif-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 

c  setup  forecast,  open  input  and  output  files 
c 

if  (debug  )  write (*,*)  'SURF' 


ydetail 
sur  f  y 
surf 
brklO 

j-ii 

ihtll 
ihtl2 
xshif t 
xdelt_gr 


'y' 

. true . 

. true . 
.false . 

1 

wspd 
wdir 
0 . 0 

2 . 0*dcal 


ccccccccccccccc 

open(41,  file  =  'dSd.list',  status  =  'old') 
do  irun  =  1,80 

read(41, ' (a) ',  end=42)  file_in 
write (*, *) file_in 
if  (irun  . gt .  0.0)  then 
ccccccccccccccc 


c 

write  ( * , * ) 

'Enter  Input  File  Name 

c 

read (*, ' (a) 

' )  file_in 

c 

file_in  =  ' 

'd3d_99.dat ' 

c 

i=0 

j=0 

do  while  (i  .eq.  0) 
j  =  j  +  l 

if (f ile_in ( j : j )  .eq.  '.')  then 
fend  =  j 
i=l 
endif 


enddo 

c 

file_out  = 

f ile_in ( 1 : fend-1 )  /  /  ' 

.  out  ' 

c 

f ile_dat 

=  file_in ( 1 : fend-1 ) // 

' dat . dat 

c 

file_in2  = 

f ile_in ( 1 : fend-1 )  /  /  ' 

par . dat ' 

CLH 

cydist=f ile_in ( 1 : fend-1 ) 

read(cydist,  ' a4 , i4 ' ) cdum,  iydist 

write  ( * , * 

)  'ydist=',  iydist 

c 

ydist=iachar (f ile_in (5 : fend-1 

)  ) 

C 

f ile_in2 

=  'D101016.dat' 

f ile_in2 

=  ' input . par ' 

c 

write ( * , * ) 

'file_in  =  ',file_in 
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c 

write ( * , * ) 

'file_out  =  ',file_out 

c 

write ( * , * ) 

'file_dat  =  ',file_dat 

c 

write ( * , * ) 

'file_in2  =  ',file_in2 

c 

iunit=20  in 

common . inc 

open ( iunit 

,  f ile=f ile_out,  status= ' unknown 

c 

write ( * , * ) 

file_out, '  OPENED' 

c 

open ( 30 , 

f ile=f ile_dat,  status= ' unknown ' ) 

c 

write ( * , * ) 

file_dat, '  OPENED' 

open (11,  f ile=f ile_in2,  status= ' old ' ) 

c 

write ( * , * ) 

file_in2, '  OPENED' 

read  (11,*)  gainina2,  per,  wspd,  wdir,  xdelt 
close  (11) 

c  write (*,*)  file_in2, '  CLOSED' 

fqd=l/per 

call  c_gainina  (gainina2 ,  igamma) 
j gamma  =  igamma 

close  (10) 

open  (10,  file  =  file_in,  status= ' old ' ) 
c  write  (*,*)  file_in,  '  OPENED' 

j_ii  =  -1 
do  j=l, 1000 

read  (10, *, end=lll)  xxin ( j ) , zzin ( j ) , tmp ( j ) , 
tmpv ( j ) , tmptheta ( j ) 

tmptheta(j)  =  tmptheta (j)  *  pi  /  180.0 
c  write (*,*) xxin ( j ), zzin ( j ), tmp ( j ) ,  tmpv ( j ), tmptheta ( j ) 

ndepth= j 
enddo 

111  continue 

close  (10) 

CLH 

xmin=xxin (ndepth) 
c  write (*,*)  ndepth 


C  re-grid  to  even  spaced 

call  c_regrid (ndepth, xxin, zzin,  xdelt_gr , nnnl , xxl ,  dxy) 

call  c_regrid (ndepth, xxin, tmp,  xdelt_gr , nnn2 , xxl ,  tmph) 

call  c_regrid (ndepth, xxin, tmptheta, xdelt_gr, nnn3, xxl, thetatemp) 
call  c_regrid (ndepth, xxin, tmpv,  xdelt_gr , nnn4 , xxl ,  v) 

c  write ( * , * ) ' nnn  =  [ ' , nnnl , ' ,  ' , nnn2 , ' ,  ' , nnn3, ' ,  ' , nnn4 , ' ] ' 

nnn  =  nnnl 

do  ii=l,nnn 

dp  =  dxy  (ii) 
c  get  wave  number 

xkd  =  0.01 

call  wavenum ( f qd, dp, xkd) 
xk  =  xkd 
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10  =  2 . 0*pi  /  xk 

hrms  =  tmph (ii) /I . 42 

c  write(*,*)  dxy(ii),  tinph(ii),  hrms,  thetatemp  ( ii )  ,  v(ii) 

xoff  =  xxl ( ii ) -xshif t 
c 

c  write(*,*)  'xoff  =  ',xoff 

c  write(*,*)  'xshift  =  ',xshift 

call  get_brk(ii,  nnn,  xxl,  dxy,  xdelt_gr,  hrms,  10, 

+  per,  xoff,  rk,  bl,  brklO,  distmax,  p) 

xtemp(ii)  =  xoff 

xktemp ( ii )  =  xk 

htemp(ii)  =  hrms 

ptemp ( ii )  =  p ( 4 ) 

iimax  =  ii 


c 

write ( * , * ) ' xoff 

= ' ,  xoff 

c 

write ( * , * ) ' xk 

X 

II 

c 

write (*,*) 'hrms 

= ' ,  hrms 

c 

write  ( * , * )  ' p ( 4 ) 

=  ',  P(4) 

c 

surf  zone  width 

index 

if  (j_ii  .it.  0 

. and .  p  ( 4 ) 

. ge .  10 ) j_ii 

enddo 

c 

write  (*,*)  ' j_ii  = 

' , j_ii, ' iimax  =  ' , iimax 

c  calculate  surf  zone  parameters,  print  to  output  file, 
c 

call  zonel(j_ii,  iimax,  dxy,  xtemp,  htemp,  ptemp,  thetatemp, 
+  xktemp,  V,  distmax,  vmax,  vmin,  thetamin,  thetamax, 

+  suml,  width,  j,  k,  hlmax,  h2max,  wid_ii) 

c  write  (*,*)  'width  =  ',widthc 

c  write (*,*) 'hlmax  =  ', hlmax 

c  write (*,*) ' h2max  =  ' , h2max 


c  15  16  18 

c  write(*,*) 'k  =  ',k,  'width  =  ', width, 'suml  =  ',  suml 

c 

c  calculate  short  output  values 
c 

call  shortcut (wdir,  wspd,  j,  iimax,  dxy,  xtemp,  suml,  k, 
ihlmax,  h2max,  per,  pet,  thetamin,  thetamax,  vmax,  vmin, 
+width,  igamma,  bl,  rk,  htemp,  wid_ii, 

ijgamma,  alfa,  bravo,  chrlie,  echo,  foxtrt,  golfl,  golf2, 
+ihtll,  ihtl2) 

c  write (*,*)  'echo  =  ',  echo 

c  write  short  output 

call  surf cast  (pet ,  depname,  Indname,  slope,  ydepth,  alfa, 
ibravo,  chrlie,  echo,  f oxtrt , golf 1 ,  golf2, 

+ihtll,  ihtl2) 
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write  ( * , * ) 

'echo  =  ' 

,  echo 

write  ( * , * ) 

'alfa  =  ' 

,  alfa 

write  ( * , * ) 

'chrlie  =  ' 

,  chrlie 

write  ( * , * ) 

'pet  =  ' 

,  pet 

write  ( * , * ) 

'echo  =  ' 

,  echo 

write  ( * , * ) 

' foxtrt  =  ' 

,  foxtrt 

write  ( * , * ) 

' jgamma  =  ' 

,  jgamma 

write  ( * , * ) 

'ihtll  =  ' 

,  ihtll 

write  ( * , * ) 

'ihtl2  =  ' 

,  ihtl2 

c  write  modified  surf  index 

CLH  add  MSI 

call  mdsrflnew  (alfa,  chrlie,  pet,  echo,  foxtrt,  jgamma, 
+ihtll,  ihtl2,  f ile_out, srfmod) 
tmpzone (irun) =golf2 
tmpxmin (irun) =xmin 
tmpydist (irun) =iydist 
tmpmsi (irun) =srfmod 
do  i=l , 4 

tmppet (irun, i) =pct (i) 
enddo 

iline=irun 

tmpfile (irun) =file_in 

write  (*,*)  'MSI  =  ',  srfmod 
close  ( 30 ) 

CLH  if  (ydetail  .eq.  'y')  call  prt_out3 (file_dat) 

close  ( iunit ) 
c 

c  end  program  surf. for 
endif 

42  continue 

enddo 
close  (41) 

open(25,  f ile= 'msi . dat ' ) 
do  ii=l,  iline 

C  write (25,  " (a, i3, 2f 7 . 1 ) " )  tmpfile(ii),  ii,  tmpmsi(ii), 

tmpzone ( ii ) 

write (25,  "(7f7.1)"  )  tmpydist (ii) , tmpmsi (ii) , tmpxmin (ii)  , 

+  tmpzone (ii) , tmppet (ii, 1) ,  tmppet (ii, 2) ,  tmppct(ii,3) 
end  do 
close  (25) 
stop 
end 
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